load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/substring.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/WRF_dyy_contributed.ncl"

; open csv file to generate the polygon
if (.not. isvar("csvfile")) then
    csvfile = "/public/users/wrf/Data/SanXia/cy_bdy.csv"
    print("csvfile is not specified, using default")
end if
print("csvfile = " + csvfile)

if (.not. isvar("region_name")) then
    region_name = "CunYi"
    print("region_name is not specified, using default")
end if
print("region_name = " + region_name)

polygon = readAsciiTable(csvfile,2,"double",1)
print(polygon)
dim_polygon = dimsizes(polygon)

if (.not. isvar("wrfout_file")) then
    wrfout_file = "/public/users/wrf/wrfout/2012111420/wrfout_d02_2012-11-14_12:00:00"
    print("wrfout_file is undefined, using default")
end if
print("wrfout_file = " + wrfout_file)

filename = systemfunc("basename " + wrfout_file)
domain = substring(filename, 7, 9)

; define output file
if (.not. isvar("output_filename")) then
    output_filename = "/public/users/wrf/Data/SanXia/cy_mask_d02.csv"
    print("output_filename is not specified, using default")
end if
print("output_filename = " + output_filename)

; open wrfout file to get the domain lat/lon information
a = addfile(wrfout_file + ".nc", "r")
xlong = wrf_user_getvar(a, "XLONG", 0)
xlat  = wrf_user_getvar(a, "XLAT", 0)
dim_domain = dimsizes(xlong)
print(dim_domain)
lonlat_recs_pool = new((/dim_domain(0)*dim_domain(1), 2/), "float")
ij_recs_pool = new((/dim_domain(0)*dim_domain(1), 2/), "integer")
rec_count = 0

; determine the max and min i,j
lon_bdry = (/min(polygon(:,0)), max(polygon(:,0))/)
lat_bdry = (/min(polygon(:,1)), max(polygon(:,1))/)
print(lon_bdry)
print(lat_bdry)

ijCorner1 = wrf_user_ll_to_ij(a, min(polygon(:,0)), min(polygon(:,1)), True) - 1 - 1
ijCorner2 = wrf_user_ll_to_ij(a, max(polygon(:,0)), max(polygon(:,1)), True) - 1 + 1

print(ijCorner1)
print(ijCorner2)

do i = ijCorner1(0), ijCorner2(0), 1
    do j = ijCorner1(1), ijCorner2(1), 1
        print("i= " + i + " j= " + j)
        if (is_point_in_polygon((/xlong(j,i), xlat(j,i)/), polygon, dim_polygon(0))) then
            print("point ("+tostring(xlong(j,i))+", "+ tostring(xlat(j,i)) + ") is in the region.")
            lonlat_recs_pool(rec_count, :) = (/xlong(j,i), xlat(j,i)/)
            ij_recs_pool(rec_count, :) = (/i, j/)
            rec_count = rec_count + 1
        end if
    end do
end do

lonlat_recs = lonlat_recs_pool(0:rec_count-1, :)
ij_recs = ij_recs_pool(0:rec_count-1, :)

output_table_body = sprintf("%7.3f", lonlat_recs(:,0)) + ", " + sprintf("%7.3f", lonlat_recs(:,1)) \
    + ", " + sprinti("%4i", ij_recs(:,0)) + ", " + sprinti("%4i", ij_recs(:,1))
output_table_all = array_append_record(region_name, output_table_body, 0)

asciiwrite(output_filename, output_table_all)
